A Chebychev propagator for inhomogeneous Schrodinger equations 
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A propagation scheme for time-dependent inhomogeneous Schrodinger equations is presented. 
Such equations occur in time dependent optimal control theory and in reactive scattering. A formal 
solution based on a polynomial expansion of the inhomogeneous term is derived. It is subjected 
to an approximation in terms of Chebychev polynomials. Different variants for the inhomogeneous 
propagator are demonstrated and applied to two examples from optimal control theory. Convergence 
behavior and numerical efficiency are analyzed. 



I. INTRODUCTION 

Inhomogeneous time-dependent Schrodinger equa- 
tions, 



ihj t \m) 



H\^{t)) + G{t)\y{t)), 



(1) 



arise in many formal solutions of quantum dynamics. In 
particular they have been employed in a time dependent 
treatment of reactive scattering [l], 0| and in optimal con- 
trol theory using time-dependent targets 0,13, 0, @] or 
state-dependent constraints 0. In reactive scattering, 
the inhomogencity results from the application of a pro- 
jection operator This projector divides the Hilbert 
space of the reactive system into subspaces correspond- 
ing, respectively, to the reactants and to the products. A 
reduced description for only the products can be derived 
where the time-dependent Schrodinger equation contains 
an inhomogencity, i.e. a source term that corresponds to 
the creation of the products 0]. 

In optimal control theory (OCT), the inhomogeneity 
may be caused by a projection operator as well. For ex- 
ample, a partitioning of the Hilbert space is implemented 
by a projection operator in order to suppress population 
in a forbidden subspace 0)- This leads to a formulation 
of OCT with a state-dependent constraint containing the 
projection operator. As a result the backward propaga- 
tion of the OCT equations includes an inhomogeneity 
in the Schrodinger equation. This term corresponds to 
the suppression of probability amplitude in the forbidden 
subspace. 

Generally, an inhomogeneous Schrodinger equation 
arises in OCT if a time-dependent target or a state- 
dependent constraint are utilized. In the common ver- 
sions of OCT, see e.g. Refs. [j| [lfj, the target is not 
explicitly time-dependent, it depends only on some final 
time T. The constraints enforce the Schrodinger equa- 
tion and a minimization of the field energy. But for ex- 
plicitly time-dependent targets 0, H, [H, @ or for a state- 



dependent constraint 0, ITll ] the optimization functional 
contains a contribution of the form 



g [ip(t),ip* (t)}dt, 



where the state of the system enters at each time t. 

For the solution of the standard homogeneous time- 
dependent Schrodinger equation, a number of numerical 
propagation schemes exist [HI, E3] ■ The Chebychev prop- 
agator [3] offers the advantage of a numerically exact 
solution. The accuracy of the calculation is then deter- 
mined by the machine precision of the computer and the 
error is uniformly distributed. The propagator is based 
on approximating the formal solution of the homogeneous 
time-dependent Schrodinger equation, 



m+dt)) = e-T^t m)h 



(2) 
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by a scries of Chebychev polynomials. Time-dependent 
inhomogeneous Schrodinger equations have been solved 
to date with split-propagator schemes 0] or via a full 
diagonalization of the Hamiltonian Q. While the lat- 
ter method is numerically expensive and quickly becomes 
unfeasible with increasing system size, the first is of only 
limited accuracy. 

Here, we derive a formal solution of the time-dependent 
inhomogeneous Schrodinger equation and we adapt it to 
the Chebychev propagation scheme. We apply this new 
propagator to the optimal control with a state-dependent 
constraint and with a time-dependent target. The paper 
is organized as follows: Section [II] presents the formal 
solution of Eq. |T]). Propagation schemes for the formal 
solution are derived in Section ITTTl The Chebychev prop- 
agation scheme is applied to OCT with a state-dependent 
constraint where the system is forced to remain in a sub- 
space of the total Hilbert space in Section IIVI In this 
case the operator G in Eq. |T]) is independent of time, 
G(t) = G. In Section [V] a second application, OCT with 
a time-dependent target, is studied keeping the full time- 
dependence, G(t). In both applications, the convergence 
of the new Chebychev propagator is analyzed in detail. 
Section [VT] concludes. 
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II. FORMAL SOLUTION 

The inhomogeneous Schrodinger equation, Eq. ([1]), is 
treated as an ordinary differential equation. It can be 
rewritten (setting h = 1) 

i||^)) = A|^)> + !$(*)>, (3) 

where |$(£)) = G(t)\(p(t)). Eq. ^ is solved subject to 
the boundary conditions 

\i>(0))=m, |$(0)) = |$o). (4) 

is known globally in the propagation time interval, 
[0, T], for example by a numerical representation on N t 
sampling points. |$(t)) is assumed to be analytic, so that 
it can be interpolated to any arbitrary point within [0, T]. 
The representation of |$(t)) on Nt sampling points corre- 
sponds to an expansion in N t basis functions. Choosing 
equidistant sampling points yields a Fourier representa- 
tion. A high-order (usually N t 1) polynomial expan- 
sion is obtained when the sampling points are chosen 
as roots of polynomials, implying non-equidistant time 
steps. The optimal representation treating correctly the 
boundaries is obtained by choosing Chcbychcv polyno- 
mials. 

The basic idea consists in devising a short-time inte- 
gration scheme for the interval [0, t] (or, more generally 
[t n ,t n+ i] with n,n+l < Nt) by taking the following 
polynomial expansion of the inhomogeneous term, 

tn — 1 m— 1 , a 

i*(*)> ^E^i^^ETfi^)- ( 5 ) 

3=0 3=0 ■>■ 

Here, Pj denotes the Chcbychcv polynomial of order j 
with expansion coefficient and f is a rescaled time, 
t = =y — 1, for t 6 [0,t]. Note that the sum on the 
right-hand side of Eq. ([5]) is only of the form of a Tay- 
lor expansion for later manipulations but the approach 
itself is not based on a Taylor expansion. The Cheby- 
chev expansion coefficients can be obtained e.g. by 



a cosine transformation (cf. Ref. [15[ and Section IIII Bl 
below). Once the coefficients are known they are 
used to generate the coefficients j n the second sum 

of Eq. ((5|) , cf. Appendix [A] This procedure has the ad- 
vantage of employing a uniform approximation of |$^^) 
in the interval [0, t}. However, if the function |4>(f)) is not 
known analytically, it needs to be interpolated at sam- 
pling points t £ [0, t] in order to calculate the expansion 
coefficients As a simpler alternative, the function 

is expanded in a Taylor series in the time interval 
[0,t]. Then the coefficients \^^) become the jth deriva- 
tive of |$(f)), at the beginning of the interval. At this 
point the properties of the global Nt interpolation func- 
tion of |$(t)) are used to calculate |$"') as numerical 
derivatives at the beginning of each interval. 

Based on Eq. ([5]), the solution of Eq. (J3j) can be written 

as 

ra-l j 

W*))(m) = E7?l A °' ) ) + ^|A^ ) }. (6) 
3=0 3 ' 

In Eq. ©, the subscript m denotes the order of the so- 
lution. The \X^') arc obtained iteratively, 

|A (0) ) = |Vo>, (7) 
|A (J) ) = -iH|A (j '- 1) ) + |$ (j '- 1) ), 
1 < j < m , 

and F TO is a function of H given by 

p m = / m (A) = ha)-™ L-™ e ^r^j 

(8) 

Equations ©-([5]) represent the formal solution to the 
inhomogeneous Schrodinger equation. 

In order to verify that Eq. © is indeed a solution to 
Eq. ([3]), let's take the derivative of Eq. ©, 



dt 



W)) 



(m) 



my- 1 



j 

m-2 



m-2 



^ * | A C*+D> + (_iA)(_iA)-™ e- iHi -E 



(-mty 



|A< 



m)\ 



3=0 



3=0 



3 



Inserting Eq. ([7]) in the first term and resumming with an upper limit m — 1, one obtains 

^ m— 1 



E ,-i 

3=0 ■ y ' 



-iH|AW) 
iA(-iA)' 



(m- 1)! 



iA|A< 



m-l)\ 



|$(™- 



-iHt 



E 

j=0 



J' 



(-iAt)j {-iHty 



(m- 1)! 



|A (r 



r 



Recognizing F m . cf . Eq. |8]) , as part of the last term and 
inserting |A' J '), cf. Eq. (J7]), in the second term, which 
then cancels with the second summand of the last term, 
this can be rewritten 



0_ 

dt 



m—1 



m)) 



(m) 



ifi £!l| A &) ) + F m |A<' 



m-l 3 

+ E> (j) >- 

3=0 ■>' 

The expression in parenthesis corresponds to \ip(t)}, cf. 
Eq. ©. Finally, replacing the second term by |$(t)), 
cf. Eq. ([S]), the inhomogeneous Schrodinger equation, 
Eq. (O, is obtained. Therefore Eq. ^ presents indeed a 
solution to Eq. ([3]). 



III. PROPAGATION SCHEMES 

The formal solution, Eq. (|6|), is subjected to a spectral 
approximation, 



\m) 



(m)N 



3=0 



(9) 



where Pm is a polynomial of order N approximating 
F m = / m (H). For example, to first order, m = 1, the 
formal solution is given by 



(i) 



-iHt 



|^o) + (-*H)- 1 (e- 



iHt 



l)|*o). (10) 



In principle one might seek a spectral approximation for 
each of the terms. However, it is numerically more effi- 
cient to rewrite the first order solution, 



(i) 



IV'o)+/i(H)(-iH|Vo) + | ( I>o; 



(11) 



with/^H) = HA)-!( 



-iHt 



1), such that only a single 



Chebychev expansion plus one extra application of the 
Hamiltonian are required. Similarly, to second and third 
order the formal solutions can be written 



m))m = \^)+t\\^) 

m))m = l^d) + *|A«) 

+/ 3 (H)|A( 3 )) 



/ 2 (H)|A (2 



(12) 



(13) 



with / 2 (H) and / 3 (A) given by Eq. © and |A ( ^) by 
Eq. ([7]). The strategy is then to seek a polynomial ap- 
proximation for the functions /,-(A). For /o(A) = e~ ?Ht , 
this corresponds to the standard Chebychev propagator 
for the homogeneous Schrodinger equation [14j . It will 
be briefly reviewed for clarity, followed by a discussion of 
the polynomial approximation of the new functions. 



A. General idea of the Chebychev propagator 

The Chebychev propagator is based on treating the for- 
mal solution as a function of an operator which is applied 
to some state vector, 

k>> = /(A)h/>>, 

and to approximate this function by an expansion in 
Chebychev polynomials P n , 

/(A)h/>) = 5>„p„(AM. 

n 

Since the Chebychev polynomials arc defined within the 
range [—1, 1], the Hamiltonian has to be renormalized, 



H, 



, H — E m i n 1 

' AE 



1 



where E m in denotes the smallest eigenvalue of H and 
AE = E max — E m i n the spectral range of H. The wave- 
function propagated from time zero to t is then obtained 
by 



N 



\m) 



i(iAE +Emin )t a n P n (-iH norm )\i>o 



n=0 



where the phase factor in front of the sum is due to the 
renormalization. 

The algorithm to estimate \<f>) = f(H)\tp) proceeds as 
follows: 



1. Calculate the expansion coefficients a n , 



2 - &n f 1 f(x)P n {x) 



dx . 



(14) 



For the function f(x) = fo(x) = e~ lxt , the integrals 
can be solved analytically resulting in the Bcsscl 
functions 14 1. 
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2. Calculate P n (— ii^norm)\ipo) using the recursion re- 
lation of the Chebychev polynomials, 



|0o> = \ipo) , 

\4>i) = -iH norm |V>o) i 

\4>n) = — 2iH raorTO |^)„_i} 



(15) 



£>n-2/ 



3. Accumulate the result, taking into account the 
phase factor due to renormalization, 



/v 



\m) 



e -i(±AE+E min )t 



n=0 



Task 1 has to be performed only once, while 2 and 3 are 
repeated for each propagation step. The number N of 
Chebychev polynomials is chosen such that the coefficient 
a^+i becomes smaller than some specified error e. Since 
the coefficients can be determined analytically, e may 
correspond to the machine precision of the computer. 



B. Chebychev propagator for inhomogeneous 
equations 

According to Eqs. (|10M13[) a single Chebychev expan- 
sion plus a few applications of the Hamiltonian are re- 
quired to solve the inhomogeneous Schrodinger equation. 
In particular, the application of a function / m (H) to a 
state vector |A ( ^ m - ) ) has to be considered, 



/™(H)|A (r 



(16) 



where |A^ m ') is obtained recursively by Eq. {7}. The 
Chebychev propagation scheme now consists in calcu- 
lating |A( m >) and performing tasks 1-3 for / m (H)|A^) 
instead of J (H) I A ( " l) ) . There are two differences with 
respect to the standard Chebychev propagator, i.e. be- 
tween the approximation of /o(H) and the approximation 

of /,„>o(H). 

(i) For m > 0, the integrals required to calculate the 
expansion coefficients a n , Eq. (|14p . cannot be solved an- 
alytically anymore. They can, however, be obtained nu- 
merically with sufficient efficiency and accuracy fl5T ] . The 
approach of Ref. [l5| is repeated here for completeness: 
Applying a Gaussian quadrature to the Chebychev poly- 
nomials, the integral in Eq. Q14[) is rewritten, 



2 -6 N_1 

= TT^— ^ fm(Xk)Pn(%k) , (17) 



N 



Since the Chebychev polynomials can be expressed in 
terms of cosines, P n (x) = cos(n9) with 9 = arccos(a;), 
Eq. (fT7|) corresponds to a cosine transformation, 



2 - 8„ 
N 



N-l 



^2 fm{6k)cos(n9 k ) . 



(19) 



fc=0 



The expansion coefficients are thus most easily evaluated 
by Fast Cosine Transformations. 

Special care is, however, required for small values of 
the argument of / m (H) since / m (H) involves division by 
the mth power of H. It is recommended to employ the 
definition of / m (H), Eq. ([8|), only if the argument is larger 
than some small value e and to use a Taylor expansion 

Of / m (H), 



/m>o(H) — f' 



f^ {]+m)\ 



(20) 



for arguments smaller than e. 

(ii) As explained in Section [TH a Taylor expansion of 
|3>(t)) for each integration interval [0, t] is most easily 
used if the inhomogeneous term is represented numeri- 
cally at sampling points. Recursive calculation of |A^ m ^) 
then requires numerical evaluation of the time derivatives 
of the 'inhomogeneous state vector', |$(t)). Second order 
m = 2 requires the first derivative which can be obtained 
with sufficient accuracy by Fast Fourier Transformation 
(FFT) and multiplication in frequency domain. However, 
an error is introduced due to finite values of | <£(£)) at the 
boundaries of the grid, i.e. t = and t = T. This error 
is increased and propagated throughout the time interval 
in the calculation of higher order derivatives by FFT and 
multiplication in frequency domain. 

Higher order schemes require therefore a different 
method for the evaluation of the derivatives. A suitable 
choice represents |$(t)) by an expansion into Chebychev 
polynomials. The derivatives are then calculated recur- 
sively based on the analytical properties of the Cheby- 
chev polynomials [l6|. Note that this implies a non- 
equidistant time grid since the roots of the Chebychev 
polynomials need to be taken as sampling points, 



T ( ( niT 

r tn = 2( 1 + cos (tvw 

n = 0,...,N t -l, 



k=0 



where the time interval t G [0, T] is scaled to t' G [— 1, 1] 
and N t denotes the number of sampling points. A non- 
equidistant time grid requires calculation and storage of 
the Chebychev expansion coefficients of the propagator 
for each At n . This additional effort is well paid off since 
a higher order scheme allows for larger time steps, i.e. 
smaller Nt. 



where the sampling points x k correspond to the TV roots 
of the Polynomial Pjy, 



C. Summary of the Algorithm 



Xk = cos 



The algorithm to solve the inhomogeneous time- 
dependent Schrodinger equation is summarized in order 
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to outline the flow chart of the numerical implementa- 
tion: 

1. Set the time grid {t n } and determine the inhomo- 
geneous term |$(£„)). 

2. Calculate the Chebychev expansion coefficients of 
F m = /m(H), cf. Eq. (H|). This needs to be done for 
each time step size that occurs in a time grid with 
non-equidistant steps or only once if an equidistant 
time grid is employed. 

3. Determine all required in Eqs. © and (J7J). 
This is done differently whether the uniform ap- 
proximation or the Taylor expansion over the short- 
time interval [0, t] is employed. 

Uniform approximation (Chebychev expansion): 

(i) Determine Chebychev sampling points {t^} 
within each short-time interval [0, t] (or, gen- 
erally [f n ,t n+ i]), cf. Eq. (fig)) , and determine 
the value of |$(rfc)) analytically or by numer- 
ical interpolation. 

(ii) Calculate the Chebychev expansion coeffi- 
cients of |$(Tfe)) by cosine transformation 
analogously to Eq. (|19[> with f m replaced by 

I*)- 

(iii) Derive for all {£„} fr om I $ , ) using 
the transformation between a Chebychev and 
a Taylor expansion given in Appendix [Aj 
i.e. employing the recursive relations given in 
Eqs. and ([AT)) . 



Taylor expansion: 

Calculate as time derivatives of |$(f n )) either 

by FFT for equidistant time grid or by numerical 
differentiation based on the analytical properties 
of Chebychev polynomials [l6[ for non equidistant 
time steps. 

4. Perform the time propagation, i.e. for each time 
step from to t (or, generally, t n to i n +i), 

(i) determine |A^) according to Eq. 0, 

(ii) obtain F m |A ( --')) by Chebychev recursion, cf. 
Eq. (T5]), 

(iii) evaluate \^{t)) (or, generally, \^>{t n+ i))) ac- 
cording to Eq. ([6|). 

In order to avoid storage of the Chebychev expansion 
coefficients, step 2 in the case of a non-equidistant time 
grid and step 3 in the case of the uniform approximation 
may also be performed during the time propagation in 
step 4. 

The algorithm to solve the inhomogeneous time- 
dependent Schrodinger equation requires adjustment of 
two parameters - the order of the algorithm, to, and the 
propagation time step, i.e. the size of the short-time in- 
terval. The latter can also be specified in terms of the 



overall number of time steps, N t , which is particularly 
useful for non-equidistant time steps. The two parame- 
ters are linked to each other: Employing a higher order 
to allows for taking larger time steps, or, equivalently, 
reducing Nt ■ A detailed discussion how the two parame- 
ters should be chosen is presented for the applications in 
Sections llV] and IV] below. 



D. 



Approximation of the formal solution by a 
Taylor expansion 



A simplified version of the algorithm is devised by ap- 
proximating the formal solution explicitly by a Taylor 
expansion. This yields a standard Chebychev propaga- 
tor plus additional terms involving derivatives of | $(£)}. 
An existing Chebychev propagation code then needs lit- 
tle modification to be adapted to an inhomogeneous 
Schrodinger equation. The accuracy of the simplified al- 
gorithm is, however, limited by the Taylor expansion. 

As shown in Appendix 151 the formal solution, Eq. (JS]), 
can be written as 

m— 1 

|^(i)) M =e- iAt |^o)+E^+i|^ (j) > ! (21) 
i=o 

with 



"j+i 



HH) 



-CB-i) 



-iHt 




Taking the Taylor expansion of the exponential in F J+1 
up to order j + 1, 



3+1 



(—i(it) k 

fci 



E 



fc=0 fc=0 

one obtains an approximation of Eq. (|2f p . 



(-iHt) k (-mt) j+1 

~ fc! - + (7 + 1)1 ' 



m-1 j+1 



j=0 



(j + 1)! 



(22) 



To all orders, only the standard Chebychev propagation 
scheme plus calculation and storage of the derivatives of 
the inhomogeneous term, |$^), are required. For exam- 
ple, the approximate solution to second order reads 



\m)(2) « /o(H)i^o) 



t |$(o)) + L|$(D 



(23) 



Unlike the propagator described in Sections IIII Bl and 
IIII C[ this scheme is explicitly based on a Taylor expan- 
sion. It is therefore expected to be valid only for small 
time steps. In that low order scheme (e.g. m = 2) 

is sufficient. The necessary derivatives can then be cal- 
culated by FFT and multiplication in frequency domain. 
The validity of this approximation is discussed in detail 
in Section [TVBj 



6 



IV. APPLICATION I: CONTROL WITH A 
STATE-DEPENDENT CONSTRAINT 

In our first example, the operator occuring in the in- 
homogeneous term is not time-dependent, G(t) = G. 

A. Model 

A simplified model of the vibrations in a Rb 2 molecule 
is considered taking into account three electronic states 
7] . The corresponding Hamiltonian reads 

H = ^H,® | ei )( ei | +£e(*)- (|ei)<e a |+ (24) 
i=i ^ 

|e 2 >(ei| + |e 2 )(e 3 | + |e 3 )<e 2 |^ , 

where Hi denotes the vibrational Hamiltonians, Hi = T+ 
Vi(R), p, the transition dipolc operator, assumed to be in- 
dependent of R, and e(t) the electric field. The electronic 
states are associated to X x I]+(5s + 5s), 1 S+(5s + 5p) and 
1 H g (5s + 4d) , with the potentials found in Ref. fi~7| . 

OCT is tested for the objective of population transfer 
from the vibrational ground state of the electronic ground 
surface to a particular vibrationally excited state via Ra- 
man transitions between the ground and the second elec- 
tronic surface. For strong laser fields e(t), such as those 
found by OCT algorithms, population at intermediate 
times will be excited not only to the second but also to 
the third electronic surface. This is particularly the case 
for the electronic states of our example, where transition 
frequencies and Franck-Condon factors for transitions be- 
tween |ei) and {e^} are very similar to those for transi- 
tions between |e 2 ) and \e3). Assuming that the third 
electronic state corresponds to a loss channel, e.g. an 
intermediate state in resonantly enhanced multi-photon 
ionization or an autoionizing state, population transfer 
into this state should be avoided at all times. This can 
be communicated to the OCT algorithm by formulating 
a statc-dcpcndcnt constraint which maximizes the pro- 
jection onto the allowed subspace 0, i.e. onto electronic 
states 1 and 2. The complete functional for optimization 
is then given by 

J[tp,tp*,e] =J q [(pt,<Pt} + Ja[e] + Jb[<p,<p*] (25) 

with 

Jo[<PT,<pi] = AoMT)|0|v>(T)>, (26) 

J a [e] = f \ a (t)[e(t)-e r (t)] 2 dt, (27) 
Jo 

Jb[f,^] = [ A^WIPsiiowlvK*)}*- ( 28 ) 

Jo 

The operator D in Jo is given by the projector onto the 
target level in the electronic ground state. The state- 
dependent constraint contains the projector onto the al- 
lowed subspace, P a iiow ln R-°f- 0j propagation via full 



diagonalization of the Hamiltonian was employed and 
only 11 vibrational levels in each electronic state were 
considered. Representing the vibrational Hamiltonians 
on a Fourier grid [l2T | and utilizing the inhomogeneous 
Chebychev propagator, the full potentials can be taken 
into account (A^id = 128). 

B. Test of the propagator 

In order to test the new Chebychev propagator, our 
numerical results are compared to those obtained by the 
symmetrical method for the Hamiltonian comprising of 
11 levels in each electronic state 0. The details of the 
symmetrical method are reviewed in Appendix [UJ The 
inhomogeneous Schrodinger equation reads 

j t \m) = -\ A[ C (t)] \m) + aAuowM*)) , (29) 

with the "initial" condition 

m = T)) = -\ D\<p(T)). (30) 

In OCT, \tp(t)) is the backward propagated wavefunction 
which is coupled to \f(t)) obtained by forward propaga- 
tion, 

j t W(t)) = - l -H[e{t)]\m (31) 

with \tp(t = 0)) = \tfo}' The convergence with respect 
to the time step At and the order m of the Chebychev 
method will be discussed in Section IIV CI 

We start by propagating the initial state, v = of 
the electronic ground state, forward with a Gaussian 
pulse given by Eq. (33) of Ref. 0. The inhomoge- 
neous Schrodinger equation is then propagated backward 
with the same field. All parameters arc taken to be 
equal to those of Ref. ■ Figure Q] shows the expecta- 
tion value of the projector onto the forbidden subspace, 
(P forbid) (t), for normalized \ip(t)). Results for the sym- 
metrical method and the first order Chebychev propa- 
gators based on Eq. ([6|) and on Eq. (|22|) (Taylor) are 
compared. A good agreement between the Chebychev 
propagators and the symmetrical scheme is found. The 
time step needed to obtain converged results with the ap- 
proximate solution is, however, smaller by a factor of ten. 
This is not surprising since the approximate solution is 
based on a Taylor expansion which requires a very small 
At. 

Our main motivation for the development of the modi- 
fied Chebychev propagator lies in its application in OCT 
with a state-dependent constraint or a time-dependent 
target. Therefore the new propagator is also tested in an 
iterative solution of the control equations. The control 
problem is that of Ref. 0. However, in the present work 
the full vibrational Hamiltonians are employed. 

Convergence of the OCT algorithm is measured by the 
normalized functional, J norm = J/(\q + ^bT), approach- 
ing one as the control equations are itcratively solved, cf. 
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FIG. 1: (color online) Normalized expectation value of the 
projector onto the forbidden subspace for a single backward 
propagation with a Gaussian pulse. Comparison between the 
symmetrical scheme of Ref. [31 with At = 1 a.u. (black solid 
line), the Chebychev propagator in first order with At = 1 a.u. 
(orange dashed line) and the Chebychev propagator in first 
order based on the Taylor approximation, Ed .([22 [I. with At = 
0.1 
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FIG. 2: (color online) J n0 rm as a function of the number of 
iterations. Solution obtained by the symmetrical scheme of 
Ref. [3] with a time step At — 1 a.u. (black solid line), the 
Chebychev propagator for first order with a time step At = 
1 a.u. (orange dashed line) and the second order approximate 
solution of Eq. (|22|l with a time step At = 0.1 a.u. (blue dotted 
line). 



Fig. [21 The symmetrical scheme (At = 1 a.u., black solid 
line) is compared to the Chebychev propagator in first 
order (orange dashed line) and to the Chebychev propa- 
gator based on the Taylor approximation, Eq. (|22|). (blue 
dotted line) in second order. An overall good agreement 
is found. The difference between the results at interme- 
diate iterations is attributed to the different models, the 
full vibrational Hamiltonian in the case of the Chebychev 
propagators and the Hamiltonian consisting of 33 levels 
in the case of the symmetrical method. 



order m 


At 


Nt 


iVcheby 


applications of H 


CPU time 


1 


2 a.u. 


180.000 


6 


1.260.000 


12 s 


2 


4 a.u. 


90.000 


8 


900.000 


lis 


3 


6 a.u. 


60.000 


10 


780.000 


9s 



TABLE I: CPU time required for one backward propagation 
to obtain converged solutions in first, second and third order 
for a total propagation time of 8 ps with equidistant time steps 
At. Also listed are the number of time steps, Nt, the number 
of times the Hamiltonian is applied and the number of terms 
in the Chebychev expansion, iVcheby 



Application of the Chebychev propagator based on the 
Taylor approximation, Eq. ([2"2")l , in first order failed: The 
inhomogeneous Schrodinger equation is solved with in- 
sufficient accuracy and the property of monotonic con- 
vergence of the Krotov method is lost during the 
iterative solution of the control equations. This is easily 
rationalized in terms of the very fast oscillations occuring 
in the optimized field. They require cither a very small 
time step or a higher order in the approximation of the 
exponential, e _lHAt , by a Taylor expansion. The com- 
plexity of the pulse and hence the time variation of the 
inhomogeneous term determine the required order of the 
Taylor approximation. 



C. Convergence behavior 

The convergence and the efficiency of the inhomoge- 
neous Chebychev propagator are analyzed with respect 
to the number of time steps Nt and the order m of the 
solution. The main numerical effort is required for the 
application of the Hamiltonian and the calculation of the 
derivatives. For a given propagation time T one would 
like to identify optimum values of N t and m that yield a 
mininum computation time for a specified accuracy. In 
general, decreasing the number of time steps N t or, re- 
spectively, increasing At will require a larger number of 
Chebychev polynomials in the expansion of the function 
f m (H ) , but also a higher order m of the solution. The re- 
cursive calculation of the Chebychev polynomials in the 
expansion of / m (H) implies continued application of the 
Hamiltonian. Moreover, higher order solutions require 
additional applications of the Hamiltonian, cf. Eqs. (fTUl 
I13[) . and determination of derivatives of the inhomoge- 
neous term up to degree m — 1 . 

We consider first equidistant time steps and evaluate 
the time derivatives by FFT. Figure [3] traces the popu- 
lation in the forbidden subspace to illustrate the conver- 
gence behavior for dynamics under a Gaussian pulse. In 
first order, cf. Fig. [Sh., the dynamics are found to be con- 
verged for time steps At < 2 a.u. At larger time steps, 
deviations are observed, in particular at short times (oc- 
curing late for backward propagation). Dynamics with 
the largest possible time step at each order are shown in 
Fig. [3b. 
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FIG. 3: (color online) Normalized expectation value for the 
projector onto the forbidden subspace demonstrating the con- 
vergence of the inhomogeneous Chebychev propagator: (a) for 
first order solutions (m = 1), convergence is lost by increasing 
At, (b) converged results for the largest possible time step at 
a given order m. 



In order to decide whether it is numerically more ef- 
ficient to keep a low order demanding a small time step 
or to employ a higher order allowing for a larger time 
step, Table U compares the CPU time required to obtain 
converged solutions in first, second and third order. The 
number of applications of H includes both those occuring 
in the Chebychev recursion and those due to the addi- 
tional terms in Eqs. (|10til3|l . For example, for third order 
and At = 6 a.u., 10 Chebychev polynomials are sufficient 
to approximate /3(H) Eq. (fT5)) . Each time step then 
implies thirteen applications of H, ten for the Cheby- 
chev recursion plus three for the additional terms (one 
for |AW) and two for |A^)), cf. Eqs. ((T3J) and ©. As 
can be seen in Table HI in terms of CPU time, it is more 
efficient to employ a higher order solution. In the con- 
text of OCT calculations, in addition to saving computa- 
tion time, a higher order propagator also allows for sav- 
ing memory since the backward propagated wavefunction 
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FIG. 4: (color online) Normalized expectation value of the 
projector onto the forbidden subspace for a Hamiltonian with 
slow time-dependence (RWA) employing a non-equidistant 
time grid and calculating time derivatives in terms of Cheby- 
chev expansions for different orders m. Results are shown for 
the smallest possible number of time steps, Nt, at given order 
m except for m = 2, Nt = 18000 (blue dashed line) which 
illustrates a non-converged case. 



needs to be stored for each time step. An inherit limit to 
increasing the time step is, however, posed by the time- 
dependence of the Hamiltonian. Expressing the formal 
solution of the homogeneous time-dependent Schrodingcr 
equation by the exponential, e~ MAt , assumes H to be 
constant within the time interval At. In our example 
the upper limit at third order, At = 6 a.u., is due to the 
breakdown of this assumption for the forward propagated 
wavefunction, \cp(t)), entering the inhomogeneous term, 

AfePallowI^C*)}- 

In order to demonstrate that a higher-order solution for 
the inhomogeneous Schrodinger equation allows indeed 
for a large time step, we modify our example by invoking 
the rotating- wave approximation (RWA) . This eliminates 
highly oscillatory terms from the field, e(t), keeping only 
the time-dependence of the envelope which is several or- 
ders of magnitude slower. However, when increasing the 
time step and the order, numerical determination of the 
derivatives by FFT and multiplication in frequency do- 
main breaks down, cf. Section [III Bl Accurate numerical 
calculation of the time derivatives of the inhomogeneous 
term, |$(t)) = XbP a ii OVJ \ip(t)) , is afforded by expanding 
|<J>(t)) in Chebychev polynomials. The sampling points 
of the time grid then need to be chosen as the roots of 
the Chebychev polynomials, leading to non-equidistant 
time steps (dividing the time interval into equidistant 
time steps corresponds to a Fourier representation) . Fig- 
ure Hand Table [n] illustrate the convergence behavior for 
a Hamiltonian with slow time-dependence. All results in 
Fig. [3] arc shown for the smallest possible number of time 
steps except the blue dashed line (to — 2, N t — 18.000) 
which illustrates a non-converged case, cf. the deviation 
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order m 


N t 




CPU time 


2 


45.000 


12.6 a.u. 


170 s 


3 


18.000 


31.4 a.u. 


69 s 


4 


9.000 


62.8 a.u. 


35 s 



TABLE II: CPU time required to obtain converged solutions 
in second, third and fourth order for a total propagation time 
of 8ps with non-equidistant time steps. Also listed are the 
smallest possible number of sampling points for the time grid, 
Nt, and the corresponding maximum time step, At max . 



v = 2 



from the converged results at short times. The number 
of time steps can be significantly reduced by employing 
higher-order schemes. The evaluation of the derivatives 
by Chebychev expansion is, however, more costly, leading 
to overall larger computation times than in Table |TJ It is 
obvious from Table [II] that this variant of the inhomogc- 
neous Chebychev propagator will unfold its full power for 
a time-independent Hamiltonian that occurs e.g. in reac- 
tive scattering calculations where a very large time step 
together with a high-order scheme will be numerically 
most efficient. The fact that the permissible time step 
can be increased by employing a higher-order solution 
illustrates that inhomogeneous Chebychev propagator is 
based on a global representation. 



V. APPLICATION II: CONTROL WITH A 
TIME-DEPENDENT TARGET 

In our second application, the operator occuring in the 
inhomogeneous term, G(t), is explicitly time-dependent. 



A. Model 

In principle it should be possible to prescribe by a laser 
pulse an arbitrary pathway that the quantum system 
should follow. To this end, OCT with a time-dependent 
target has to be employed [J, In the total func- 
tional, Eq. (|25[) . the final-time term then disappears, 
Jo[tpT, < / 5 tX = 0, and the state- and time-dependent term 
becomes Q 

J b [<p tV fi] = I \ b {<p(t)\G(t)\<p(t))dt. (32) 
Jo 

Maximization of J b corresponds to X b > and fulfills the 
conditions for monotonic convergence [tJ. 

A simple model comprising of five of the 33 levels of 
Section IIVI are taken to mimic a double A-system, cf. 
Fig. [5] Initially all population is assumed to be in v = 0, 
and at the final time the population in v = 2 is to be 
maximal. Additionally, the time interval [0, T] is divided 
into subintervals where the population of the interme- 
diate levels v' = 6, v = 1, and v' = 7, is maximized, 
i.e. we prescribe a 'trajectory' where the ladder of the 
double A-systcm is sequentially climbed up. While this 



v=0 



v= 1 



FIG. 5: Prescribed 'trajectory' for a time-dependent target: 
climbing up the ladder of a double A-system. 



represents a simple toy model, it serves the purpose of 
illustrating the case where the operator of the inhomoge- 
neous term of the Schrodingcr equation, G(t), is explicitly 
time-dependent. The inhomogeneous equation for back- 
ward propagation reads, 

j t \m) = -\ Aft*)] \m) + A 6 e(t)b(t)> , (33) 



with the "initial" condition 



W = T)) = 0. 



(34) 



Dividing the time interval [0, T] into four subintervals, 
< Ti < T 2 < T 3 < T, the target is defined as the 
projector onto v' = 6 in [0, Ti], onto v = 1 in [Ti,T 2 ], 
onto v' = 7 in [T 2 ,T 3 ] and onto v = 2 in the subinterval 
[7s, T], 

G(t) = |6)(6|9(Ti-t) + 

\l)(l\Q(t-T 1 )0(T 2 -t) + 

\7)(7\e(t-T 2 )0(T 3 -t) 

|2)(2|0(t -T 3 )0(T-t) (35) 

with O(i) the Heaviside function. In order to avoid nu- 
merical problems due to discontinuities, &(t) is approxi- 
mated by 



i 



-ki 



(36) 



where the parameter k determines the steepness with 
which the target level changes. While for large values 
of k, the step function is recovered, small values of k 
imply overlap in time of two different targets near the 
Ti, i = 1 — 3. In the following k is varied between 
k = 1CP 4 a.u. and k = 10 4 a.u. The final time is set 
to T = 5.4 ps. The subintervals are taken to be of the 
same length, Ti = 1.35 ps, Ti = 2T X and T 3 = 3Ti. 
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FIG. 6: (color online) Time-dependent target: Time evolution 
of the level populations with an optimal field (shown in grey). 
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11130 


11072 


11014 
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11172 


11114 


11057 



TABLE III: Transition frequencies in cm 1 for the five levels 
employed as time-dependent target. 



B. Results 
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The new propagator is tested for the iterative solution 
of the control equations with the time-dependent target. 
The guess field consists of a sequence of four 7r-pulscs, 
one in each time interval. Figure [6] shows the evolu- 
tion of the level populations using the optimized field for 
k = 10 4 . They follow by and large indeed the prescribed 
'trajectory'. Population of levels other than the target 
one and fast oscillations in the populations are observed 
only when switching from one target to the next. The 
spectrum of the optimized field is shown in Fig. [7] The 
transition frequencies of our model, listed in Table IIIII 
and indicated in Fig. [TJd, arc contained within the spec- 
trum. Additional frequencies which do not correspond 
to the main transition frequencies are observed. They 
are attributed to the complexity of the optimal solution 
which may include beatings between levels, Stark shifts 
etc. 

The improvement of the time-dependent target func- 
tional with the number of iterations is demonstrated in 
Fig. [5] for different values of the steepness parameter k. 
Monotonic convergence is observed. However, the algo- 
rithm cannot reach 100%. We attribute this to the way 
the target is switched and overlap in time of different tar- 
gets is created around the J*, i = 1 — 3. For large k the 
changes in the target functional are almost instantaneous 
and cannot be followed by the dynamics, cf. the oscilla- 
tions in the level populations in Fig. [6j However, at the 
same time, the targets do almost not overlap. This yields 
the highest value of the target functional, about 85%. 
For smaller values of k the dynamics can follow more 



FIG. 7: (color online) (a) Spectrum of the optimized field, 
(b) Detail of the spectrum (w £ [10975 cm~\ 11205 cnT 1 ]). 
The numbers and arrows indicate the six main transition fre- 
quencies of the model, cf. Table IITT1 



smoothly. However, the overlap between different tar- 
gets is increased, i.e. contradictory objectives are asked 
at the same time. This decreases the value of the target 
functional to about 79%. 



C. Convergence behavior 

The convergence of the inhomogeneous Chebychev 
propagator is again analyzed with respect to the time 
step At and to the order m. We restrict ourselves here 
to the case of equidistant time steps and calculation of 
the derivatives by FFT and multiplication in frequency 
domain. The convergence behavior is illustrated in Fig. [5] 
by the time evolution of the final target level (v = 2) pop- 
ulation. Converged results are obtained for At < 4 a.u. 
(At < 10 a.u.) in first (second) order, i.e. a larger time 
step than in Section IIVI can be used. We attribute this 
to the much simpler model. 

Table [TV] compares the CPU time required to obtain 
converged solutions in first and second order. The same 
conclusion is obtained as in Section IIVI i.c it is more 
efficient to employ a higher order scheme. 

Overall, no difference in the convergence behavior for 
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FIG. 8: (color online) The renormalized functional Jt, cf. 
Eq. (|32p . as a function of the number of iterations for pa- 
rameters k corresponding to different overlaps in time of two 
targets. 
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FIG. 9: (color online) (a) Time-dependent population of the 
final target level v = 2: Converged results in first and sec- 
ond orders for the largest possible time step (black solid and 
orange dotted lines). Also shown is a non-converged result in 
first order (blue dashed line) - the deviations from the con- 
verged solutions are evident in the insert. 



time-dependent and time-independent operators in the 
inhomogeneous term, G(t) is found. This can be ratio- 
nalized as follows: the convergence is determined by the 
fastest timcscale of the dynamics, i.e. by the rapid os- 
cillations of the field. The time-dependence of the pro- 
jection operator introduces a time-dependence which is 
much slower and hence does not affect convergence. 



order m 


time step At 


Afchcby 


applications of H 


CPU time 


1 


4 a.u. 


8 


495.000 


0.72 s 


2 


10 a.u. 


10 


264.000 


0.62 s 



TABLE IV: CPU time required to obtain converged solutions 
in first and second order for a total propagation time of 5.4 ps 



VI. CONCLUSIONS 

A formal solution to the time-dependent inhomoge- 
neous Schrodinger equation was derived based on an 
expansion of the inhomogeneous term. Three levels of 
Chebychev approximations are involved. 

(i) The first one yields the Chebychev propagator 
where the argument of the Chebychev polynomi- 
als is the Hamiltonian. Truncating the expansion 
at the desired order m, the formal solution is sub- 
jected to a spectral representation with Chebychev 
polynomials. A propagation scheme similar to the 
standard Chebychev propagator for homogeneous 
Schrodinger equations is then obtained: Instead of 
e -iHAt a f unc tion / m (H) is expanded in Chebychev 
polynomials. For the exponential function, the ex- 
pansion coefficients can be calculated analytically, 
for the function / m they need to be obtained numer- 
ically. This is achieved by Fast Cosine Transforma- 
tions, utilizing the definition of Chebychev polyno- 
mials in term of cosines. 

(ii) The second level expands the inhomogeneous state 
vector |<&(£)) in Chebychev polynomials within each 
short-time integration interval [0, t\. The argument 
of this expansion is the rescaled time t covering the 
interval. This Chebychev approximation is easily 
applied only if is known analytically. If |$(t)) 
is determined numerically on sampling points cov- 
ering the global propagation time interval [0,T], 
there are two choices. | <&(£)) needs to be interpo- 
lated to sampling points within [0, £]. In a simpler 
alternative, the Chebychev expansion is replaced 
by a Taylor expansion based on numerical deriva- 
tives at the beginning of each time step. This has 
been done for the present applications. 

(iii) The numerical calculation of the derivatives re- 
quires a third level of Chebychev approximation 
where the argument is the time t covering the global 
propagation time interval [0, T]. This implies a 
non-equidistant time grid where the derivatives are 
evaluated according to the procedure described in 
Ref. [l6|]. This expansion overcomes the numerical 
error introduced by non-zero boundary values of 
the inhomogeneous state vector |<E>(t)} at £ = and 
t = T . An alternative based on equidistant time 
steps employs Fast Fourier Transforms and multi- 
plication in frequency domain. However, in that 
case the errors introduced at the boundary of the 
time grid build up. Therefore this scheme is limited 
to low order where only first or second derivatives 
are required. 

An even more approximate solution to the time- 
dependent inhomogeneous Schrodinger equation is ob- 
tained by rewriting the formal solution explicitly in terms 
of a Taylor expansion. The propagator then consists of 
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the standard exponential term plus time derivatives of 
the inhomogeneous terms. This approximation is numer- 
ically less efficient than the propagator for the full for- 
mal solution. Moreover, it may become instable in opti- 
mal control applications where the inhomogeneous term 
often is highly oscillatory and the numerical evaluation 
of derivatives by FFT becomes difficult. The main ad- 
vantage of this propagation scheme lies in the fact that 
it requires very little modification of existing standard 
Chebychev propagation codes. 

Both Chebychev propagation schemes were tested in 
two optimal control applications. OCT with a state- 
dependent constraint [7(, e.g. maximizing population in 
an allowed subspace of the Hilbert space, yields a time- 
independent operator in the inhomogeneous term while 
an explicitly time-dependent operator is obtained in OCT 
with a time-dependent target [1, 0, H, @ . Convergence 
of the propagation schemes was demonstrated for both 
applications. The convergence behavior was studied in 
detail as a function of the order of the solution and the 
required number of time steps for a given overall propaga- 
tion time. For applications with a fast time-dependence 
such as OCT, a low order scheme with a small time step 
and evaluation of the time derivatives by FFT was found 
to be the best choice. For applications with a slow time- 
dependence or time-independent Hamiltonians such as 
reactive scattering calculations where large time steps are 
permissible, a high-order scheme is numerically most effi- 
cient. This reflects that the propagation scheme is based 
on a global representation of the inhomogeneous term. It 
is this regime where the new propagator can best unfold 
its power. 

The new Chebychev propagator provides a stable and 
accurate numerical solution to the time-dependent inho- 
mogeneous Schrodinger equation. It is most efficient for 
high order and large time steps. Ideally an inherent time- 
dependence of the Hamiltonian should also be incorpo- 
rated into the Chebychev scheme. This is the subject of 
a further study. 
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APPENDIX A: TRANSFORMATION TO OBTAIN 
|$ (J) ) FROM 

When the solution of the inhomogeneous Schrodinger 
equation is based on the uniform approximation, a trans- 



formation linking the Chebychev expansion coefficients, 
to the coefficients of the Taylor expansion, \$^}, 
cf. Eq. ([5]), is required. In other words, given a 
vector [Ao, . . . A m ] T we want to compute the vector 
[B 0l . . . B m ] T such that 

m rn j 

Y,A m ,kPk(x) =Y,B m A , (Al) 

k=0 3=0 

i.e. we identify and \& 3 ^) to A rn ^ and B m j respec- 
tively. Let 

P fc (aO = E C ^4- ( A2 ) 

3=0 J ' 

Since Chebychev polynomials obey the recursion relation, 
P k+l {x) = 2xP k {x) - P k -i{x) , (A3) 

we obtain 

fc+l j k k — l j 

X^-'A » -y- (A4) 

3=0 J ' 3=0 J ' 3=0 J ' 

or 

k+1 x 3 k+1 x 3 kl X 3 

r ''— = 2 J2 Ck >3-i (j _ lv -J2 Ck -h3^- 

3=0 J ' 3=1 KJ >' J=0 ■'■ 

(A5) 

Hence, the C coefficients satisfy 
Ck+ifi = — Cfc-i,o , 

Cfe+ij = 2jCk,j-i - Cfc-i,3 , 1 < j < k - 1 

Cfc+i.fc+1 = 2(fc + l)C fe , fe . (A6) 

Based on this result we can compute the B coefficients 
recursively, 

Bi+ij = Bij + A i+ i ti+ id + ij , l<j<i 
-Si+i.i+i = Ai y iCi y i , (A7) 

for 1 < i < m — 1 with 

Bo,o = A),o j Bi,o = ^i,o j B ly i = . (A8) 

APPENDIX B: PROOF OF THE EQUIVALENCE 
OF EQS. © AND (\21\) 

It is shown by induction that the formal solution, 
Eq. ©, and Eq. (|2"Tj) are equivalent. Writing Eqs. © 
and (|2ip for m = 1, one obviously obtains in both cases 
the equation of the first order, Eq. (flO]) . We assume that 
Eqs. © and (j2~lj) arc equivalent in order m — 1, i.e, 
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2 P 



j=o J J=0 



\m)(m-i) = E -l AW ) + ^-ilA^) = e- m \i> ) + E F 3 - +1 |*fc>) , (Bl) 

3=0 J ' 

let us now prove that they are equivalent in order m. 

ro-i j 

W*)>(m) = E 7T|A (i) ) + F m |AW) 



7 

3=0 J 

m-2 ,,• ,m-l 



3=0 

We continue by showing that F rn _i = F m (— «H) + 1 



y-3 fr>L—±. , 

E -yl^) + 7^— Tv l A(m_1) ) + ^ (hh)^™- 1 )) + |$(» 

E -y|A^>> + FJ^™- 1 )) + f F m (-iH) + 1 (^Yyy J lA^) (B2) 



(m-l)! ' 
m— 2 



F m (-iH) + 1 = HH)-<—) e"** - E ^ - + iT^T 

(m-l)! \ p< Q j\ (m-l)! I (m-l)! 



3=0 
m-2 



3=0 7 



m—2 ,a 



= (- l H)-(™- 1 ) ( E I - F m _r . 

Equation (|B2[) thus becomes 

M*)>(m) = E L -\ Xi ^ + Fm-llA*™- 1 ') + F m |$(— !)) 
3=0 J ' 

Making use of our assumption. Eq. (jBip . we obtain 

\m)(m) = e-^'^o) + E p J+il* W) ) +F r „i$ (rn - l) ) 

3=0 
m — 1 

= e-^l^+E^'+il^)- 

3=0 

This concludes the proof. 



m-2 



APPENDIX C: ALGORITHM OF THE Schrodinger equation, Eq. ©, 

SYMMETRICAL METHOD 



The symmetrical method of Ref. is based on the m) = e -iAt ^ (0)) + g -iH* / ^ =. (r))dr _ (C1) 



formal integration of the inhomogeneous time-dependent 
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Assuming |$(r)) to be constant in [0, t] and taking its 
value to be 

|t(T)) _ mm±mi weM , ,C2) 

the integral in Eq. (|C1[) can easily be computed and we 
obtain 

m)) « e- iA *|^(0)) + HH)- 1 (e-' lAt - l) |$(0)> . 

(C3) 

This solution is formally equivalent to the first order of 
the Chebychev propagator, cf. Eq (flO]). However, the 



evaluation of \ip(t)) proceeds differently in the symmet- 
rical method and the first order Chebychev propagator. 
The latter subjects Eq. (|C3[) to a spectral approxima- 
tion. It only requires a representation of the Hamilto- 
nian such that its action on a state vector can be evalu- 
ated. A numerically very efficient representation is based 
on the Fourier grid where evaluation of H|V>) scales as 
0(N log N) with N the number of grid points. The sym- 
metrical method diagonalizes H(£) at each time step in 
order to directly employ Eq. (|C3[) . Since diagonalization 
scales as 0(N 3 ), where iV is the dimension of the Hilbcrt 
space, this is feasible only for sufficiently small N. 
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